#ifndef _REGENERATION_ENGINE_HPP_
#define _REGENERATION_ENGINE_HPP_

#include "definitions.h"
#include "engine.hpp"

#include <fstream>
#include <iostream>
using namespace std;

#include <cmath>
#include <boost/math/distributions/students_t.hpp>
using namespace boost::math;

#define ALPHA 0.05
#define MIN_REGC_LEN 20
#define REG_CYCLES_MIN 30


class Regeneration_Engine {
public:
  Regeneration_Engine(Engine &e)
  : engine(e)
#ifdef LOG
  , logfile(*new ofstream("reg_averages.dat", ios_base::app))
#endif
  {}

  ~Regeneration_Engine() {
#ifdef LOG
    delete &logfile;
#endif
  }

  void init(Queue *dq, Queue *nq, unsigned long nj, unsigned long m) {
    reg_dq = dq;
    reg_nq = nq;
    reg_nj = nj;
    reg_mcl = m;
  }

  void run_pilota(Queue *cycle_queue, vector<Job> &jobs) {
    bool reg_cond = false;
    double A = 0;
    unsigned long v = 0;
    Job *moving_job;
    Queue *departure_queue, *arrival_queue;
    double event_time;
    do {
      //fire 1 event
      engine.step(departure_queue, arrival_queue, moving_job, event_time);
      if(arrival_queue == cycle_queue) {
  	//update A and v
  	++v;
  	A += moving_job->completion(event_time);
      }
      reg_cond = check(departure_queue, v, jobs);
    }
    while(!reg_cond);
    cout << "length of the reg. cycle: " << v << endl;
    cout << "point estimation: " << A / v << endl;
  }

  void run(Queue *cycle_queue, double thickness_target, vector<Job> &jobs) {
#ifdef LOG
    cc = vector<unsigned long>(jobs.size(), 0);
#endif
    n_reg_cycles = REG_CYCLES_MIN;
    unsigned long n_reg_cycles_old = 0;
    //accumulators
    double A=0, Sv=0, Saa=0, Svv=0, Sav=0;
    //run
    bool go_on = true;
    do {
      unsigned long nc = n_reg_cycles - n_reg_cycles_old;
      for(unsigned long i=0; i<nc; ++i) {
#ifdef LOG
	cerr << "> starting reg. cycle #" << n_reg_cycles_old + i << endl;
#endif
	bool reg_cond = false;
	double Ac = 0, vc = 0, event_time;
	do {
	  //fire 1 event
	  Queue *departure_queue, *arrival_queue;
	  Job *moving_job;
	  engine.step(departure_queue, arrival_queue, moving_job, event_time);
	  if(arrival_queue == cycle_queue) {
	    //update A and v
	    ++vc;
	    Ac += moving_job->completion(event_time);
	  }
	  reg_cond = check(departure_queue, vc, jobs);
	}
	while(!reg_cond);
#ifdef LOG
	cerr << "completions (ccj): ( ";
	for(unsigned int j=0; j<jobs.size(); ++j) {
	  cerr << jobs[j].get_n_completions() << " ";
	  jobs[j].reset_completions();
	}
	cerr << ")" << endl
	     << vc << " completions\n---" << endl;
#endif
	//update accumulators
	A += Ac;
	Saa += Ac * Ac;
	Svv += vc * vc;
	Sav += Ac * vc;
	Sv += vc;
#ifdef LOG
	logfile << Ac << endl;
#endif
      }
      //compute INTERVAL
      point = A / Sv;
      double tq = quantile(students_t(n_reg_cycles - 1), 1 - ALPHA/2);
#ifdef LOG
      cerr << "t(" << n_reg_cycles-1 <<"," << 1 - ALPHA/2 << ") = " << tq << endl
	   << "r: " << point << " Saa: " << Saa << ", Sav: " << Sav << ", Svv: " << Svv << endl;
#endif
      offset =
	tq *
	sqrt((double)n_reg_cycles / (n_reg_cycles - 1)) *
	sqrt(Saa - point * Sav * 2 + point * point * Svv) / Sv;
      //compute THICKNESS
      thickness = offset * 2 / point;
#ifdef LOG
      cerr << "interval: " << point << " +- " << offset
	   << ", thickness: " << thickness*100 << " %"
	   << " / " << thickness_target*100 << " %" << endl
	   << "job cycles per reg. cycle (avg): " << (double)Sv / n_reg_cycles << endl;
#endif
      //check halt condition
      go_on = thickness > thickness_target;
      if(go_on) {
	//more reg. cycles
	n_reg_cycles_old = n_reg_cycles;
	n_reg_cycles = (unsigned long)ceil(pow(thickness / thickness_target, 2) * n_reg_cycles_old);
#ifdef LOG
	cerr << "new number of reg. cycles: " << n_reg_cycles << endl;
#endif
      }
    }
    while(go_on);
    events_per_cycle = (double)Sv / n_reg_cycles;
#ifdef LOG
    cout << "n. reg. cycles: " << n_reg_cycles << endl
	 << "job cycles per reg. cycle (avg): " << events_per_cycle << endl;
#endif
  }

  void get_interval(interval_t *res) {
    res->confidence = 1.0 - ALPHA;
    res->point = this->point;
    res->offset = this->offset;
    res->thickness = this->thickness;
    res->n_reg_cycles = n_reg_cycles;
    res->events_per_cycle = events_per_cycle;
  }

private:
  Engine &engine;
  double point, offset, thickness;
  unsigned long n_reg_cycles;
  double events_per_cycle;
#ifdef LOG
  vector<unsigned long> cc;
  ofstream &logfile;
#endif

  //reg. cond. control variables
  Queue *reg_dq;
  Queue *reg_nq;
  unsigned long reg_nj;
  unsigned long reg_mcl;

  bool check(Queue *departure_queue, unsigned long v, vector<Job> &jobs) {
    bool reg =
      (reg_dq == departure_queue) && //departure from terms (useless?)
      (reg_nq->get_n_jobs() == reg_nj); //(N, 0, 0)

    bool len = (v >= reg_mcl); //reached min. length of a reg. cycle

    bool jobs_cond = true;
    if(reg)
      for(unsigned int i=0; i<jobs.size(); ++i)
	jobs_cond &= (jobs[i].get_n_completions() > 0);

#ifdef LOG
    if(reg)
      cerr << "REG. POINT; length: " << v << "; " << (jobs_cond ? "ok":"no") << " jobs_cond" << endl;
#endif

    return reg && len && jobs_cond;
  }
};
#endif
